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A general method for the direct evaluation of the temperature dependence of the quantum- 
mechanical reaction rate constant in many-dimensional systems is described. The method is based 
on the quantum instanton approximation for the rate constant, thermodynamic integration with 
respect to the inverse temperature, and the path integral Monte Carlo evaluation. It can describe 
deviations from the Arrhenius law due to the coupling of rotations and vibrations, zero-point energy, 
tunneling, corner-cutting, and other nuclear quantum effects. The method is tested on the Eckart 
barrier and the full-dimensional H + H2 — >■ H2-I- H reaction. In the temperature range from 300 K to 
1500 K, the error of the present method remains within 13% despite the very large deviations from the 
Arrhenius law. The direct approach makes the calculations much more efficient, and the efficiency is 
increased even further (by up to two orders of magnitude in the studied reactions) by using optimal 
estimators for reactant and transition state thermal energies. Which of the estimators is optimal, 
however, depends on the system and the strength of constraint in a constrained simulation. 



I. INTRODUCTION 

The measurement of the temperature dependence of 
the rate constant is one of the important tools of chem- 
ical kinetics in determining mechanisms of chemical 
reactionsjir— Significant deviations from a simple expo- 
nential behavior can be evidence of tunneling and of other 
nuclear quantum effects^i^ These effects are particularly 
strong for hydrogen transfer reactions with a high acti- 
vation barrier or at low temperatures. Recently, how- 
ever, quantum effects have been observed also in many 
enzymatic reactions at physiological temperatures.®^^ It 
therefore becomes more and more important to have ac- 
curate theoretical methods for computing the tempera- 
ture dependence of the rate constant. 

Probably the oldest yet still the best known expression 
for the thermal rate constant k(T) at temperature T is 
the empirical Arrhenius lawj^^ 

fc^(T) = Ae-^"/'^^'^. (1.1) 

Here is the Boltzmann constant, Ea the activation 
energy, and the temperature dependence is purely ex- 
ponential. An improvement over the Arrhenius law was 
provided by the transition state theory (TST);^^— in 
which 

where h is the Planck's constant, Q^{T) and Qr{T) are 
the partition functions of the transition state and the re- 
actants, respectively, and IS.E'^ is the barrier height for 
the reaction. Here the temperature dependence includes 
a fractional power T" in addition to the exponential. 
Nevertheless, both Arrhenius law and TST are basically 
purely classical, so they cannot take into account tunnel- 
ing and other nuclear quantum effects. 



Although the simplest quantum effects can be taken 
into account within the TST in an ad hoc fashion, 
by replacing the partition functions by their quantum 
analogs for the simple harmonic oscillator (which takes 
into account the zero-point energy, the Wigner tunnel- 
ing correction, and approximate quantization of the vi- 
brational motion), a more systematic approach requires 
quantum treatment of the nuclear motion. This is of 
course extremely difficult, and therefore various approx- 
imate yet accurate methods have been developed. These 
include, e.g., the semiclassical methods ^^"^^ or the so- 
called quantum transition state theoriesj^"— 

In this paper, we evaluate the temperature dependence 
of the rate constant starting from the quantum instan- 
ton (QI) approximation!^ This quantum transition state 
theory has been shown to describe correctly not only all 
of the above-mentioned quantum effects, but also corner- 
cutting, coupling of vibrational and rotational motions, 
multiple tunneling paths, etc. As we evaluate the tem- 
perature dependence of the rate constant directly, we can 
speed up the QI calculation significantly by avoiding the 
tedious umbrella sampling necessary for computing the 
rate constant itself. Furthermore, if it is only the tem- 
perature dependence of the rate constant that is needed, 
we can increase the accuracy of the QI approximation 
for the rate constant by canceling some small remaining 
systematic errors. This can be useful, e.g., if we know 
fc(To) at a temperature Tq very accurately and would 
like to know k{T) at other temperatures. The temper- 
ature dependence of the rate constant is computed via 
a thermodynamic integration^^''^" with respect to the in- 
verse temperature. A similar thermodynamic integration 
in the framework of the QI model was used by Ceotto 
and Miller to compute the rate constant for several one 
and two-dimensional systems using a discrete variable 
representation!^ Below, we develop a general thermody- 
namic integration procedure based on the path integral 
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implementation, which is suitable for many-dimensional 
systems. 

The method is tested on two simple systems for 
which exact quantum dynamical calculations are feasi- 
ble; the Eckart barrier and the full nine-dimensional H 
+ H2 ^ H2+ H reaction. While the former system is 
the simplest one-dimensional model of a bimolecular re- 
action, the latter is the simplest bimolecular chemical 
reaction with an energy barrier. As such, it has been 
widely studied and attained the status of a benchmark 
reactioni^^i^ Despite the apparent simplicity, this reac- 
tion remains a challenging test for new approximations. 
This is due to the presence of strong quantum effects, 
as the lightest atoms are involved in both bond breaking 
and bond formation. This reaction was investigated not 
only for the temperature dependence of its rate constant, 
but also for the kinetic isotope effect^^""— the presence 
of the geometric phase cffectf^li^ etc. 

The remainder of this paper is organized as follows: 
Section II describes the methodology, i.e., the QI approx- 
imation, the thermodynamic integration with respect to 
the inverse temperature, the path integral formalism, and 
the relevant estimators. Computational details and, in 
particular, the analysis of various errors are presented 
in Sec. III. Section IV contains the results. First, it is 
shown how the low and high temperature limits are ob- 
tained. The numerical results for the Eckart barrier and 
the H + H2 reaction are then presented. The results are 
compared with the Arrhenius law, TST, TST with the 
Wigner tunneling correction, and the exact quantum cal- 
culation. Finally, we discuss how the efficiency is related 
to the dependence of the statistical error on the number 
of imaginary time slices in the path integral, and how 
this error, in turn, depends on the system under study. 
Section V concludes the paper. 



II. METHODOLOGY 

A. Quantum instanton approximation for the 
thermal rate constant 

The QI approximation for thermal rate constants was 
introduced in Ref. i28|. The most direct derivatio n^^'^^i*^" 
starts from the exact Miller-Schwartz-Tromp formula for 
the rate constant4i 



applied to Eq. (|2.ip yields the QI approximation for the 
rate constant, '^^-'i^ 

fc(T)«.MT) = i;C.(/3,0)f ^ (2.3) 

where AH{f3) is a specific type of energy variance^ 
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The delta-delta correlation function used above is defined 
as 

Cdd W,t) = Tr (e-^^/2A,e-^^/2e^^*/'"'A,e-^^*/'"' 

(2.5) 

where the generalized delta function operator is given by 
= A K (f ) ~^^]^S K (f) - C7] m-'/^ II veil (2.6) 

and e(r) is the reaction coordinate such that ^(r-l-) = at 
the transition state. Similarly, ^(r) = defines the posi- 
tion of the dividing surface 7. We have used mass-scaled 
coordinates in which all degrees of freedom have the same 
mass TO. In practice, the exact delta function constraint 
is approximated by a Gaussian constraint corresponding 
to a harmonic constraint potential T4onstr(r)r^'^ 
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(2.8) 



The accuracy of the QI approximation has been al- 
ready verified in numerous applications.— "^i^"^"— The 
main shortcoming of the QI method is the neglect of 
recrossing which is however, neglected in any quantum 
or classical transition state theories. The recrossing ef- 
fects on the quantum instanton rate constant have been 
quantified for several collinear reactions by Ceotto and 
Milleri.2i Fortunately, the recrossing effects become gen- 
erally less important in higher dimensions. 



B. Temperature dependence via the 
thermodynamic integration 



k{T)Qr^ / dtCa{l3,t) (2.1) 
Jo 

where Cff (/3, t) is the symmetrized fiux-fiux correlation 
function, 



The goal of this paper is to compute the temper- 
ature dependence of the rate constant, i.e., the ratio 
k{T2)/k{Ti). Within the QI approximation, this ratio 
is given by 



(2.2) 

with Hamiltonian operator H , inverse temperature (3 := 
l/ksT, time t, and Fj the flux operator through the 
dividing surface 7. A stationary-phase approximation 



HT2) ^ QrWi) Ag(/3i) Cdd W2,0) Sfe°o) 
fc(Ti) AH{h) Cdd (A, 0) ' 



(2.9) 



where we multiplied and divided the numerator and de- 
nominator by Cdd(/5,0). In this expression, quantities 
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AH{f3) and Cs {(3, 0) /Cm 0) can be computed di- 
rectly by the Metropolis Monte-Carlo procedure because 
they are thermodynamic averages. On the other hand, 
the ratios Qr^/Qr^ and Cdd (/32, 0) /Qd (/3i, 0) 
cannot be computed this way since they involve ratios 
of quantities at different temperatures. These ratios can, 
however, be calculated by the method of thermodynamic 
integratio n^^ ^^'^ with respect to the inverse temperature 



potential $ is given by 
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(2.10) 
(2.11) 



where Er and E^ are the thermal energies of the reac- 
tants and of the transition state, respectively. A sim- 
ilar thermodynamic integration was used within a dis- 
crete variable representation of the QI approximation to 
compute the rate constant for several collinear triatomic 
reactions.— Unlike Qr and Cdd, the energies are normal- 
ized quantities because they can be written as logarith- 
mic derivatives: 



EriP) - 
E\I3) := - 



dlogQrW) _ dQril3)/dl3 



dl3 

dlogCdd(/3,0) 
dp 



Qr{f3) ' 

rfgdd(/3,0)/rf/3 
Cdd(/3,0) 



(2.12) 
(2.13) 



Hence they can be computed directly by a Monte Carlo 
procedure. 



C. Path integral representation of relevant 
quantities 

Quantum thermodynamic effects can be treated rigor- 
ously using the imaginary time path integral (PI).**^"^^ 
Let D be the number of degrees of freedom {D — 1 for 
the Eckart barrier and £) = 9 for the H3 potential) and 
P the number of imaginary time slices in the PI. The 
PI representations of the partition function^"— and the 
delta-delta correlation function^ are 



{f3)^C dr(i) • • • / dr(^) 



exp 



-/3$({r(^)}) 

(2.14) 

C^d 0) = C 1 dr(i) ... 1 dr(^) A (r^")) 
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where C = [mP/{2T:h'^ (3)]^-^^^ , r^*^ is a D-dimensional 
vector representing the sth time slice, and the effective 
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(2.16) 

For P = 1, the above expressions reproduce classical 
statistical mechanics, while exact quantum statistics is 
reached in the limit P ^ 00. 

In practice, there are two main strategies for evaluating 
thermodynamic averages using the PI: the PI molecular 
dynamics (PIMD>^>^ or PI Monte Carlo (PIMC);^ We 
use the PIMC procedure together with the Metropolis 
algorithm. The basic idea is to sample the PI configura- 
tion space according to an appropriate weight p, which 
is, e.g., for Cdd given by 



P 



:= A 



.(0) 



A 



(2.17) 

and then, at each sampled configuration, to evaluate the 
so-called estimator {^A^^}) of the relevant physical 
quantity A. The final estimate of A is given by the aver- 
age (^^) along the PIMC trajectory. 

Using the PI representation of Qr and Cdd, one can ob- 
tain estimators for all quantities needed in Eq. (12. 9p . i.e., 
AH, Cff/Cdd, and the logarithmic derivatives of Q^-, Cdd 
(i.e., the energies Er, E^). Those for AH and Cs/Cm 
are listed in Ref. Isil. 



D. Estimators for Er 

The simplest estimator for the energy Er , the so-called 
Barker or thermodynamic estimator (TE),^^ can be de- 
rived directly from Eq. (|2.12l) and the PI expression 
(HH, giving 



P DP mP^ 



(^) 



(2.18) 

As observed by Herman et al.,^^ the TE can have a large 
statistical error, which can be avoided with the so-called 
virial (VE)^ or centroid virial (CVE)^ estimators. In- 
voking the virial theorem, the kinetic energy in these 
two estimators is replaced by an expression involving the 
gradient of the potential energy This is convenient in 
the PIMD implementations since the gradient is already 
available. In PIMC simulations, however, only the po- 
tential is needed for the random walk, and in order to 
avoid computing the gradients, alternative approaches 
have been proposed. One can, e.g., employ the centroid 
thermodynamic estimator^^ or more generally, use a pro- 
cedure based on rescaling coordinates^SiS in which the 
gradients of the potential are replaced by a single deriva- 
tive that can be evaluated by finite difference.— Variants 
of the latter approach have been applied successfully to 
compute thermal energies and heat capacities, ®° kinetic 
isotope effects^i^, equilibrium isotope effects^^, or the 
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derivatives of the flux-flux correlation function^S needed 
in the generalized QI model. "^^ 

The VE for Er can be derived most directly by the 
change of coordinates x^"*) := in the PI (|2.14p . 

yielding 



-i 

s—1 



Vir 



dV[il + qy/\^'^ 



dq 



, (2.19) 



where g is a small dimensionless parameter and the q- 
derivative is evaluated by finite difference at g = 0. Sim- 
ilarly, the CVE can be obtained by the change of vari- 
ables x^*' := /3~i/2(r(*) — r^'-^^), where one first subtracts 
the so-called centroid coordinate r^*^^ := P^^ X^f^i 
The resulting estimator is 
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r.CVE 
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(2.20) 



1 ^ f^(^(.)) ^ dV[ri^) + [l + qf/\vi^)-v(^))] \ 
^ s=i dq \ 



E. Estimators for 

In the case of constrained simulations near the tran- 
sition state, the constrained weight function p.l7p can 
be approximated by using the Gaussian approximation 
of the delta function from Eqs. (I2.7p - (|2.8p . Besides a 
prefactor, this amounts to adding a constraint potential 



$constr({r^'^}) := V;onstr(r 



(P/2) 



riP) 



to the effective potential $. Assuming that V^constr is inde- 
pendent of temperature and following a derivation similar 
to that for estimators of Er, one obtains the thermody- 
namic, virial, and centroid virial estimators for E^ , 



^TE ~ ^r.TE n 'i'constr 
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(2.21) 



d$constr[(l + g)i/'{r(^)}] 



dq 



(2.22) 
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^+$co„str({r(^)}) 



CVE — ^r,CYE ^ 

d^constrKrC^) + (l-f g)i/2(r(-) - r(C))}] 



dq 



(2.23) 



Although the above estimators converge to the exact 
results, we found that the statistical errors can be de- 
creased slightly by employing an alternative set of esti- 
mators, derived using an exact relation 



,({r(^)}))=r\ 



which is valid for a harmonic constraint potential for any 
value of P. The new estimators are 



-^TE ~ ^r,TE' 

pt,P _ pP 
^VE ~ ^r,VE 



rf$constr[(l + g)^/Mr(-)}] 

dq 



(2.24) 
(2.25) 
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d$co„str[{r(^) + (1 + (7)^/'(r(^) - r(^))}] 



CVE — ^r,CVE 



dq 



(2.26) 



Estimators (|2.24p - (|2.26p are in a way more intuitive than 
estimators (|2.2ip - (|2.23p : in the limit of a sharp con- 
straint, the constrained energy should be independent 
of the type of constraint. 

It should be stressed that the last terms in the VE 
and CVE in Eqs. (g^Sl), (1^^ . (1^^ . and ^^TI^ are 
important; without them the agreement among the TE, 

VE, and CVE is lost. In other words, an intuitive guess 
t p p 

such as i?cvE ~ CVE would not give a correct answer 
for the constrained energy. 

Finally, we also tested a constraint potential that is 
proportional to temperature, i.e., 



14c 



(2.27) 



where Vconstr is a harmonic potential independent of tem- 
perature. As a result, the constraint (|2.7I) itself is ac- 
tually independent of temperature. Following again a 
derivation similar to that for estimators of one ob- 
tains another set of the TE, VE, and CVE for E^ , that 
look exactly like the estimators (|2.24l) - (|2.26l) for Konstr 
independent of temperature. The only difference is that 
the random walk is done with a different constraint. 



III. COMPUTATIONAL DETAILS AND ERROR 
ANALYSIS 

All calculations were performed with a PIMC code im- 
plemented in Fortran 90. Sampling of the configurational 
space in the PIMC simulation was done using three types 
of moves. Staging algorithm^ was employed to move 
all unconstrained beads. Constrained beads, i.e., beads 
s — P/2 and s = P which feel the constraint potential 
Konstr, were sampled with the free particle single slice 
algorithmi^ Finally, whole chain moves^ were used to 
speed up sampling of the potential energy surface. 

The Gaussian constraint potential must be strong 
enough in order to exert the constraining effect on the 
system. When this condition was satisfied, the converged 
results were independent of the constraint. However, the 
statistical root mean square error (RMSE) of the tran- 
sition state energy E^ increases with the strength of the 
constraint because sampling of the configuration space 
becomes more difficult. Therefore the selected strength 
of the constraint should take into account these two ef- 
fects. We have used A: = 10 a.u. in both systems. 
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All quantities needed in the ratio (I2.9p were evalu- 
ated using the above mentioned estimators. The ther- 
modynamic integrations (|2.10p and (|2.1ip were evalu- 
ated with the Simpson rule using 25 values of (3 between 
/3o = i/ksTo and /3max = l/fcs7inin with the reference 
temperature Tq — 1500 K and the minimum tempera- 
ture Tmin = 200 K. The number of beads was chosen to 
be inversely proportional to the temperature, with the 
maximal number of beads (used for Tmin = 200 K) being 
P = 96 for the Eckart barrier and P = 160 for the H + H2 
reaction. 

The error of the final result consists of four main error 
contributions: a) the statistical error due to the Monte 
Carlo simulation, b) the error due to the discretization 
of the TI, c) the error due to the discretization of the 
PI (i.e., the "finite P error"), and d) the actual error 
of the QI approximation. We have carefully separated 
these four contributions and attempted to make the first 
three contributions small in comparison with the error of 
the QI. In more complicated systems, this may not be 
possible and especially the final statistical error may be 
comparable to or larger than the error of the QI. Because 
the exponentiation of the TI is quite sensitive to various 
errors, a detailed analysis of errors was carried out for the 
ratio A:(200K)/fc(1500K), i.e., over the largest tempera- 
ture range, where the first three types of errors are the 
greatest. The TI was evaluated by four different numeri- 
cal methods, namely the trapezoidal, Simpson, Simpson 
3/8, and Boole methods.-S^ 

Comparing the analytical bounds on the discretization 
errors of the TI integrals using a numerical estimate of 
a higher order derivative, one can conclude that both 
the Simpson and Simpson 3/8 methods were much better 
than the trapezoidal rule and that the Boole method did 
not provide any further improvement. Specifically, for 
the Eckart barrier the error of the ratio due to the dis- 
cretization of the TI was 2%, 0.03%, 0.02%, or 0.04% for 
the trapezoidal, Simpson, Simpson 3/8, or Boole meth- 
ods, respectively. For the H -I- II2 reaction, the discretiza- 
tion error of the final ratio was 7%, 0.3%, 0.1%, or 0.3%, 
in the same order. It should be emphasized that these 
error estimates are very conservative, as the actual differ- 
ence between the final ratios based on different methods 
was an order of magnitude smaller than what one would 
expect from the error estimates. The final results dis- 
played in the plots used the Simpson method. 

The statistical RMSEs were estimated with the block 
averaging method using a variable block size^ to remove 
correlation of the PIMC data. The statistical error of the 
TI was evaluated using an appropriate formula for each 
integration method and assuming that the statistical er- 
rors of energies at different temperatures were uncorre- 
lated. As expected, the statistical error did not depend 
much on the integration method, and was always close 
to the statistical error for the Simpson method. The sta- 
tistical error of the final ratio was 0.3% for the Eckart 
barrier and 1.6% for the H-I-H2 reaction. 

The finite P error, i.e., the error due to the discretiza- 



tion of the Feynman PI, was obtained by repeating cal- 
culations of all quantities at all temperatures with twice 
smaller numbers of beads (P — >■ P/2) and then extrap- 
olating each quantity to P — >■ 00, assuming 1/P^ con- 
vergence. We emphasize that we used the extrapolated 
results only for estimating the finite P error of the com- 
puted ratio and not for estimating the ratio itself, which 
could be dangerous. The finite P error of the ratio was 
—0.3% for the Eckart barrier and —3.5% for the H + II2 
reaction. 

[We note that for H + II2 one of the temperatures 
(972.93 K, a temperature in the vicinity of which a sharp 
bend in the — dependence occurs) required a five 
times longer simulation to reduce the TI discretization 
errors. This was because a small statistical error had a 
huge effect on the estimate of the fourth derivative and 
hence on the analytical estimate of the discretization er- 
ror.] 

To sum up, in both systems the TI discretization error 
was negligible to the statistical and finite P errors, which, 
in turn, were small in comparison to the error of the QI 
approximation. 



IV. RESULTS 

A. Temperature dependence according to the 
Arrhenius law, TST, and the TST with the Wigner 
tunneling correction 

At high temperatures, the rate constant is expected 
to behave classically and follow the Arrhenius law or the 
more accurate TST result. Whereas the Arrhenius law 
(jl.ip predicts the rate constant ratio to be a simple ex- 
ponential function of the inverse temperature. 



-B„(/32-/3i) 
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TST (|1.2[) gives the ratio of rate constants 



fcTST(/32) P1QKP2) 



A:tst(/3i) P2QHPl)Qr{h) 



(4.1) 



(4.2) 



In particular, assuming the partition functions and Qr 
to be separable into products of classical rotational and 
vibrational partition functions, the temperature depen- 
dence (14. 2|) of TST rate constant includes an additional 
fractional power law besides the exponential dependence 
in the Arrhenius law (14. 1|) . 

At somewhat lower temperatures, when quantum ef- 
fects start to play a role, the basic TST expression (|4.2p 
can be improved in several ways: First, classical parti- 
tion functions Qr and can be replaced by their ex- 
act quantum analogs for a harmonic potential. Second, 
quantum tunneling can be included approximately via 
the Wigner tunneling correction. This method corrects 
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the rate constant with a multiphcative factor 



K = 1 



24 



(4.3) 



where is the imaginary frequency of the asymmet- 
ric stretch along the reaction coordinate. The correction 
can be derived by treating the motion through the tran- 
sition state as a vibration on an upside down potential 
and expanding the quantum partition function to sec- 
ond order in /3. Although an improvement over TST, 
the Wigner tunneling correction cannot describe multi- 
dimensional tunneling. 



B. Eckart barrier 

A simple model of an activated chemical reaction is 
provided by the Eckart barrier, a one-dimensional system 
described by the potential 



Vix) = Vo [cosh(aa;)]" 



(4.4) 



We use standard parameter values Vo = 1-56 • 10~^ a.u., 
a = 1.36 a.u., mass m = 1060 a.u. and reaction coor- 
dinate ^ :— X. The exact quantum rate constant /cqm 
for this reaction can be obtained by integrating the ex- 
act quantum mechanical cumulative reaction probability, 
which is known analyticallyi^ 

Figure [T] (a) compares the QI results with the exact 
QM results, TST (which is equal to the Arrhenius law 
here), and the TST including the Wigner tunneling cor- 
rection. The reference temperature is 1500 K and the 
plot shows ratios for temperatures down to 200 K. Since 
classical recrossing does not occur for the Eckart barrier, 
all TSTs should converge to the correct quantum results 
at high temperatures. The figure confirms that this is in- 
deed the case: note that all curves are tangent at the high 
temperature limit. At low temperatures, one reaches the 
quantum regime where tunneling is important and con- 
sequently the Arrhenius plot of the exact QM result has 
a large curvature. While TST has a huge error, the QI 
approximation agrees very well with the QM result. Note 
that the Wigner tunneling correction improves over the 
TST and captures the tunneling effect partially but still 
fails to recover the curvature of the exact result. 

Figure [T] (b) shows the relative error of the rate con- 
stant ratio for the different methods. Whereas both TST 
and TST with the Wigner tunneling correction deterio- 
rate rapidly with decreasing temperature, the QI method 
has an error below 3% for all temperatures above 330 K. 
The QI approximation has a significant error (> 10%) 
only at very low temperatures, below ~ 270 K. However, 
this error was well understood already in the original pa- 
per by Miller et al.— and can be remedied by considering 
two separate dividing surfaces at very low temperatures. 
(Here we have used a single dividing surface at all tem- 
peratures for simplicity.) 
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FIG. 1: Eckart barrier, (a) Temperature dependence of the 
rate constant, (b) Temperature dependence of the relative 
error of the ratio fc(r)/fc(1500 K). 



The temperature dependence of the reactant and tran- 
sition state energies is shown in Fig. [21 While both curves 
are quite smooth, small discretization errors in the inte- 
grals can have large effects on the exponentiated result. 
By a detailed error analysis described in Sec. HI, we 
found that the Simpson method was sufficient for the TI 
over the whole temperature range. Note that the VE 
for Er gives zero, but can be easily corrected with an 
analytical correction 1/(2/3). 

The three different estimators for the constrained en- 
ergy £'■!■ at T = 515.15 K are compared in Figs. [3] and 
m Panel (a) of Fig. [31 which uses the simpler estimators 
(11211)- (112SI), shows that the TE, VE, and CVE agree for 
all examined values of P and, in particular, converge to 
the same value for P — > cxd. The three estimators, how- 
ever, differ in their statistical convergence. Unlike for the 
unconstrained result, where the CVE is the optimal esti- 
mator, for the constrained energy, the optimal estimator 
is the VE. This can be clearly seen in Fig. [31 (b) which 
shows the RMSEs of the different estimators for differ- 
ent values of P. While the RMSE of the TE increases 
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with P, the RMSEs of the VE and CVE remain approxi- 
mately constant as a function of P, with the VE having a 
much smaUer statistical error. Assuming that the desired 
convergence is achieved for P = 24, the speedup factor 
achieved by using the VE compared to the TE and CVE 
is approximately 2.9^ « 8 and 8.1^ « 60, respectively. 
It is clear from the figure that both the speedup factor 
and the best estimator depend on P and hence on the 
temperature. 

Figure 2] shows the same results, but computed with 
the estimators (|2.2ip - (|2.23p . The statistical errors are 
very similar, although for the VE slightly larger than 
those in Fig. [31 



The H + Ha 



H2 + H reaction 



The temperature dependence of the rate constant 
of the H + H2 — >■ II2 -I- H reaction was studied on the 
Boothroyd-Keogh-Martin-Peterson (BKMP2) reactive 
potential energy surface<^"— The classical transition 
state of this system has a collinear configuration with 
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FIG. 3: Dependence of the transition state energy (a) and 
of the statistical RMSE of the transition state energy (b) on 
the number of imaginary time slices for the Eckart barrier 
at T = 515.15 K. The constraint potential is independent of 
temperature and estimators (|2.24|l - (|2.26|l are used. 



equal bond lengths dHaHb = '^HbHc- A suitable reaction 
coordinate is therefore given by the difference of the bond 
lengths, 



C(r) 



IHkH, 



(4.5) 



Figure [5] (a) shows the temperature dependence of the 
rate constant in the range from 200 K to 1500 K. The ex- 
act QM results are from Ref. i35j. At high temperatures 
the TST curve is tangent to the exact QM curve, but at 
low temperatures, there is a significant discrepancy even 
for the TST with the Wigner tunneling correction. On 
the other hand, the QI approximation agrees very well 
with the exact QM result all the way to 200 K. The rela- 
tive error of the rate constant ratio is shown in Fig. [5](b) 
which confirms that the error of the QI approach is within 
13% in the full temperature range whereas all other ap- 
proximations have huge errors already for temperatures 
as high as 500 K. 

In case of the H + II2 II2 + H reaction, the VE had 
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to be corrected for both Er and calculations. The rea- 
son is that the virial theorem only holds for bound sys- 
tems. The transition state of the H3 system can translate 
freely as a whole and the three translational degrees of 
freedom yield a correction of D/(2/3) = 3/(2/3) to the VE. 
In the reactant region, both the H atom and H2 molecule 
can move freely and the six translational degrees of free- 
dom give a correction of 6/(2/3) to the VE. 

The temperature dependence of the reactant and tran- 
sition state energies is shown in Fig.[S] While both curves 
are quite smooth, small discretization errors in the inte- 
grals can have large effects on the exponentiated result. 
By an error analysis described in Sec. Ill, we found that 
the Simpson method was sufRcient for the TI over the 
whole temperature range. 

Figures [7] and [8] show how the constraint energy and 
the RMSE of E^ depend on P for T = 515.15 K. Panel 
(a) of Fig. El which uses the estimators (g^ll)- (11211) 
and a constraint potential independent of temperature, 
shows again that the TE, corrected VE, and CVE give 
approximately the same results for all values of P and, 
within a statistical error, converge to the same limiting 
value for P 00. Panel (b) of Fig. 13 shows that while 
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FIG. 5: The H+ H2 -s> H2+ H reaction, (a) Temperature de- 
pendence of the rate constant, (b) Temperature dependence 
of the relative error of the ratio fc(r)/fc(1500 K). 



the statistical error of the CVE is approximately constant 
as a function of P, the RMSEs of the TE and VE grow 
with P. However, in this case, the results are quite well 
converged for P = 32 and at this point the RMSE of the 
CVE is still larger than the RMSE of the TE, although it 
is already smaller than the RMSE of the VE. While for 
lower temperatures where larger values of P are needed, 
CVE would eventually become the optimal estimator, it 
is not so for T = 515.15 K. The growth of the RMSE 
of the VE with P is due to the fact that unlike for the 
Eckart barrier, the transition state of the H3 system can 
move freely as a whole. 

Finally, Fig. |S1 shows analogous results, still using esti- 
mators (I2.24p - (|2.26l) . but obtained with a constraint po- 
tential (|2.27l) proportional to T (chosen such that the two 
types of constraints coincide for T = 515.15 K). As ex- 
pected, the statistical errors of the TE, VE, and CVE are 
similar to those in Fig. |71 [obtained with the constraint 
potential (|2.8p independent of T] . 
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A general method for the direct evaluation of the tem- 
perature dependence of the quantum rate constant was 
presented. The main advantage of this method is the 
increased efficiency: Evaluating the temperature depen- 
dence directly, without computing the rate constant at 
any given temperature, allows us to avoid a tedious um- 
brella sampling procedure. 

Besides efficiency, the direct calculation of the temper- 
ature dependence of the rate constant can also improve 
the accuracy: Our ratios fcQi(T)/fcQi(1500K) for both 
the Eckart barrier and the H -I- H2 — > H2 -t- H reaction 
have somewhat smaller relative errors than the errors ob- 
tained for the absolute QI rate constants in previous stud- 
ies of these systems j^i^ The smaller relative error in the 
ratio of rate constants is due to a favorable cancellation 
of various systematic errors, such as the systematic error 
of about 25% of the QI model at high temperatures (that 
can also be removed by an ad hoc correction of AH)2& 
and small recrossing effects in the H + H2 H2 + H re- 
action, also at high temperatures. 



It is noteworthy that for both reactions, the RMSEs 
of transition state energies depend on the strength of the 
constraint. Weakening the constraint facilitates sampling 
of the configuration space and the error of the CVE de- 
creases, approaching the well-known unconstrained sit- 
uation where the CVE is typically the optimal estima- 
tor. However, at the same time the constraint must be 
strong enough to exert the constraining effect and de- 
scribe the situation near the transition state properly. As 
a result, "ranking" of the estimators is not universal but 
can change with the potential used as a constraint and 
is in general different from the ranking for unconstrained 
simulations. 

The dependence of the error of the VE on P is best 
understood in terms of the ring polymer interpretatiort^L 
of the discretized PI: The quantum thermodynamics of 
the original system can be interpreted as the classical 
thermodynamics of the ring polymer. The constrained 
PI simulation for the one-dimensional Eckart barrier is 
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completely bound, resulting in the RMSE independent of 
P. In the full-dimensional hydrogen exchange reaction, 
on the other hand, even the constrained simulation allows 
the system to move as a whole. This is exactly where the 
VE is known to have a RMSE increasing with P . 

The CVE estimator is usually the optimal estimator 
in unconstrained systems with some translational (i.e., 
free-particle) degrees of freedom. In a system in which 
only two slices are bound (in our case, slices P/2 and P) , 
the symmetry between different slices is lost and so is, to 
some extent, the advantage of subtracting the centroid. 
This explains why the RMSE of the VE can sometimes 
be smaller than the RMSE of the CVE for all values of 
P, which we observed in Figs. [3]and|4l 

To sum up, while in generic systems at very low tem- 
peratures, the CVE is expected to be the optimal es- 
timator for energy, at finite temperatures in constrained 
simulations, the VE or even the TE can have the smallest 
RMSE. The results obtained in this paper can serve as a 
guide for choosing the best estimator for a given system. 
However, since the additional cost of evaluating all three 
estimators is negligible in comparison to the cost of the 
PIMC random walk or PIMD simulation, we recommend 
computing all three estimators, evaluating their RMSEs, 
and using the one with the smallest RMSE in a given 
situation. 
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